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ABSTRACT 


Generalized  primary/secondary  flow  equations  and  a  spatial-marching 

'  •  -  ■  \ 

solution  algorithm  have  been>utilized-to  develop  a  procedure  to  compute  the 
three-dimensional  viscous  flow  around  a  submerged  body  in  maneuver.  The  primary/ 
secondary  flow  equations  are  an  approximation  to  the  Navier-Stokes  equations  for 
flows  in  which  a  primary  flow  direction  can  be  identified.  Important  elements  of 
the  approximation  are  a  locally  specified  primary  flow  direction  and  a 
decomposition  of  the  secondary  velocity  field  to  identify  a  small  velocity  vector 
for  approximation.  No  approximations  are  introduced  for  pressure  in  this 
approach.  The  primary/secondary  flow  equations  are  a  well-posed  initial-value 
problem  in  a  spatial  coordinate  nominally  aligned  with  the  primary  flow  direction 
and  are  solved  by  a  sequentially  decoupled  implicit  algorithm.  The  procedure 
provides  an  order  to  two  orders-of-magnitude  run  time  advantage  over  solution  of 
the  Navier-Stokes  equations.  Results  are  presented  for  the  flow  past  an 
unappended  submarine  hull  in  drift  at  a  Reynolds  number  of  16  million  and 
incidence  of  20  degrees.  These  results  are  consistent  with  experimental 
observations  and  provide  a  means  to  compute  the  complex  three-dimensional  viscous 
flow  field  economically.  '  •  , 

l 
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1.  INTRODUCTION 


The  flow  around  submerged  bodies  at  different  attitudes  evokes 
considerable  interest  in  the  field  of  hydrodynamics.  Such  flow  fields  occur 
in  many  and  varied  applications.  Of  particular  interest  is  the  flow  field 
around  a  submarine,  torpedo,  or  similar  body  in  maneuver.  The  computation  of 
this  flow  field  is  important  to  the  analysis  of  performance  and  noise 
characteristics  of  the  device. 

The  flow  field  around  a  submerged  body  in  maneuver  is,  in  general, 
three-dimensional,  viscous,  and  unsteady.  The  large  scale  unsteadiness  in  the 
flow  field  is  primarily  due  to  the  time-dependent  motion  of  the  maneuvering 
body.  The  computation  of  this  unsteady,  three-dimensional,  viscous  flow  field 
with  available  computer  resources  is  a  formidable  task.  However,  many  portions 
of  a  maneuver  can  be  steady  (e.g.  a  steady  descent  or  climb)  or  quasi-steady, 
and  the  flow  field  around  the  body  can  be  obtained  from  solution  of  the  steady 
flow  equations  for  a  fixed  body  attitude.  Several  different  attitudes  can  be 
studied  in  this  manner.  The  results  would  provide  body  surface  pressures  which 
can  then  be  used  to  derive  the  forces  and  moments  on  the  body.  The  computation 
of  this  quasi-steady  flow  field  is  addressed  in  this  report. 

The  complexity  of  the  flow  around  a  submerged  body  in  steady  maneuver 
presents  considerable  difficulties  in  the  utilization  of  computational  methods 
to  compute  such  flows.  These  difficulties  arise  from  the  three-dimensionality 
of  the  flow  field  characterized  by  large  secondary  vorticity  and  velocity 
generated  that  result  in  a  flow  field  with  diverse  length  scale  flow 
structures.  Solution  of  the  three-dimensional  averaged  Navier-Stokes  equations 
avoids  making  physical  approximations  other  than  those  associated  with 
turbulence  modeling.  However,  this  approach  is  very  costly,  even  with  modern 
supercomputers,  because  the  accurate  resolution  of  the  complex  flow  structures 
requires  very  large  mesh  densities.  In  light  of  this  high  cost  of  solution, 
approximations  have  been  developed  which  reduce  the  steady  subsonic 
Navier-Stokes  equations  to  a  form  which  is  well-posed  as  an  initial-value 
problem  in  a  suitable  direction  and  could  be  solved  by  a  spatial-marching 
algorithm.  The  advantage  of  such  an  approach  is  that  spatial-marching 
algorithms  can  be  devised  which  are  considerably  less  expensive  in  terms  of 
computer  resources  (run  time  and  storage)  than  algorithms  for  the  elliptic 
Navier-Stokes  equations.  The  trade-offs  for  this  advantage  are  that  errors  are 


introduced  due  to  the  approximations  and  the  range  of  flow  problems  that  can  be 
addressed  is  restricted  relative  to  the  Navier-Stok.es  equations  because  of 
factors  such  as  flow  separation,  stagnation  points,  and  transonic  flow 
effects.  Nevertheless,  this  approach  seems  well  suited  for  a  number  of  flows 
arising  in  practical  applications  and  can  provide  very  high  resolution  of 
three-dimensional  viscous  flow  structures  at  relatively  low  cost.  In  addition, 
the  spatial-marching  approach  could  provide  a  large  number  of  detailed  flow 
calculations  at  moderate  cost  for  use  in  design  optimization  studies. 

The  potential  gain  in  economy  of  solution  for  a  class  of  complex 
three-dimensional  viscous  flows  by  utilizing  a  spatial-marching  approach 
provided  ample  motivation  to  study  and  establish  the  applicability  of  the 
approach  to  compute  the  flow  field  around  a  submerged  body  in  maneuver.  This 
study  of  the  applicability  of  the  spatial-marching  approach  to  compute  the  flow 
field  around  a  submerged  body  was  the  primary  focus  of  the  present  contract. 
Three  aspects  were  identified  and  comprised  the  study: 

(a)  Establishment  of  a  framework  of  approximations  and  governing 
equations  suitable  for  the  three-dimensional  viscous  flow  around  a 
submerged  body  in  maneuver.  The  guiding  principle  in  this  process  was 
that  the  approximations  be  valid  considering  the  physics  of  the  flow  field 
and  the  governing  equations  be  a  well-posed  initial  value  problem 
considering  mathematical  requirements  for  successful  application  of  a 
spatial-marching  solution  algorithm. 

(b)  Formulation  of  governing  equations  for  efficient  numerical  solution. 
Identification  and  application  of  efficient  numerical  solution  procedures 
to  achieve  the  run  time  advantages  desired. 

(c)  Computation  of  the  three-dimensional  viscous  flow  around  a  submerged 
body  to  demonstrate  successful  completion  of  (a)  and  (b)  and  establish  the 
applicability  of  the  approach  to  compute  the  submerged  body  flow  field. 

All  three  aspects  were  successfully  addressed  under  the  present  contract. 

The  next  section  of  this  report  describes  the  generalized 
primary/secondary  flow  equations  utilized  as  the  approximate  flow  equations  for 
the  submerged  body  flow  field  and  the  sequentially  decoupled  implicit  algorithm 
utilized  to  solve  the  equations  efficiently.  Application  of  the  procedure  to 
compute  flow  around  a  submerged  body  in  maneuver  is  also  described.  The  third 
section  of  the  report  describes  results  of  the  computations  of  the  flow  around 
a  submerged  body  in  drift  and  compares  the  results  with  available  experimental 
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data. 


The  report  concludes  with  a  description  of  current  capability  and 
suggestions  for  future  development  of  the  capability. 

2.  PHYSICAL  APPROXIMATIONS  AND  GOVERNING  EQUATIONS 

2.1  Overview 

Two  basic  types  of  physical/mathematical  approximations  have  been 
suggested  to  reduce  the  Navier-Stokes  equations  to  a  form  suitable  for  spatial 
marching  solution.  Approximations  in  both  viscous  and  inviscid  terms  in  the 
Navier-Stokes  equations  are  necessary  to  obtain  non-elliptic  (well-posed)  ap¬ 
proximate  equations.  The  viscous  approximation  entails  neglecting  terms  rep¬ 
resenting  streamwise  diffusion.  This  approximation  requires  identifying  a 
direction  in  the  flow  field  along  which  the  approximation  is  made  and  associat 
ing  the  streamwise  (marching)  coordinate  with  this  direction.  Two  types  of  in 
viscid  approximations  have  been  suggested:  (a)  an  assumed  form  for  the  steam- 
wise  pressure  gradient  term,  and  (b)  a  small  scalar  potential  approximation  fo 
the  secondary  flow.  Either  of  these  inviscid  approximations,  along  with  the 
viscous  approximation,  produces  non-elliptic  governing  equations. 

The  inviscid  approximation  which  assumes  a  given  form  for  the  streamwise 
pressure  gradient  terra  has  obvious  roots  in  two-dimensional  boundary  layer 
theory,  and  has  been  used  extensively.  The  second  type  of  inviscid 
approximation  (termed  the  small  scalar-potential  approximation)  has  been 
investigated  recently  by  Briley  and  McDonald  [1].  This  approximation  does  not 
employ  any  approximation  for  pressure  gradient  terms  and  instead  approximates 
convective  terms  in  the  secondary  flow  momentum  equation  by  neglecting  the 
scalar-potential  component  of  a  vector-decomposed  secondary-flow  velocity 
field.  The  effect  of  the  generated  streamwise  vorticity  (secondary  flow)  on 
the  pressure  field  is  retained  in  the  streamwise  pressure  gradient  with  the 
small  scalar-potential  approximation.  This  effect  is  usually  absent  in  the 
streamwise  pressure  gradient  when  it  is  assumed  to  be  of  given  form. 

The  inviscid  and  viscous  approximations  just  described  should  not  be 
confused  with  other  equations  under  the  general  heading  of  'parabolized 
Navier-Stokes'  equations  which  are  not  well-posed  initial-value  problems  but 
are  solved  iteratively  by  approximate  well-posed  spatial  marching  iterations. 
The  iterations,  when  converged,  remove  the  approximation  made  in  each  of  the 
iterates  to  obtain  solutions  of  the  elliptic  set  of  approximate  flow  equations 
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Preliminary  computations  of  the  steady  flow  around  a  submerged  body  in 
drift  with  a  spatial  marching  procedure  were  completed  by  Govindan  and  Levy 
[2],  These  computations  approximated  streamwise  pressure  gradients  in  the 
governing  equations  by  setting  them  to  zero.  At  this  level  of  approximation, 
the  results  demonstrated  that  the  qualitative  features  of  the  three-dimensional 
viscous  flow  around  a  submerged  body  could  be  computed  within  the  context  of  a 
spatial  marching  algorithm.  These  results  provided  a  basis  to  further  pursue 
and  refine  the  approach  to  provide  quantitatively  accurate  computations  of  the 
submerged  body  flow  field,  the  objective  of  the  present  contract. 

The  first  task  addressed,  under  the  present  effort,  was  to  develop  the 
framework  of  approximations  and  governing  equations  suitable  for  the  submerged 
body  flow  field.  Govindan,  Briley,  and  McDonald  [3,4]  have  recently  derived 
generalized  priraary/secondary  flow  equations  as  an  approximate  set  of  flow 
equations  applicable  to  a  wide  class  of  three-dimensional  viscous  flows.  These 
equations  are  an  approximation  of  the  steady  Navier-Stokes  equations  based  on 
knowledge  of  a  primary  flow  direction  which  includes  the  effect  of  surface 
geometry.  The  primary/secondary  flow  equations  along  with  the  small  scalar 
potential  approximation  are  a  well-posed  initial  value  problem  suitable  for 
solution  by  a  spatial  marching  algorithm.  These  equations  were  deemed  to  be 
the  preferable  set  of  equations  for  computing  the  flow  around  a  submerged  body 
in  manuever  for  two  reasons.  First,  Govindan,  Briley  and  McDonald  [4]  in  their 
computations  have  demonstrated  the  quantitative  accuracy  that  could  be  obtained 
from  the  primary/secondary  flow  equations.  Second,  the  choice  of  a  potential 
flow  model  to  generate  streamwise  pressure  gradients  for  the  pressure 
approximation  is  unclear  and,  therefore,  the  approximation  is  difficult  to 
use.  The  difficulty  does  not  arise  with  the  small  scalar  potential 
approximation  wherein  no  approximations  are  made  for  the  streamwise  pressure 
gradient. 

2.2  Primary/secondary  flow  equations  for  the  submerged  body  flow  field 

The  physical  approximations  made  in  the  primary/secondary  flow  equations 
are  a  distinguishing  feature  of  the  present  approach.  The  present 
approximations  are  similar  in  spirit  to  those  made  in  two-dimensional  slender 
channel  theory,  and  to  a  lesser  degree  in  two-  and  three-dimensional  boundary 
layers.  In  these  latter  approaches,  one  velocity  component  (the  surface-  or 
coordinate-normal  component)  is  assumed  to  be  small.  In  applications  such 
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as  flow  about  a  submerged  body,  however,  large  steamwise  vorticity  is 
generated.  The  velocity  components  associated  with  this  strearawise  vorticity 
are  known  to  be  large  (of  the  order  of  the  streamwise  velocity)  and  cannot  be 
neglected.  In  the  present  approach,  a  small  velocity  vector  is  identified 
which  provides  a  basis  for  approximating  flows  with  large  streamwise  vorticity 
and  secondary  velocity. 

Two  features  of  the  primary/secondary  flow  equations  central  to 
identifying  a  small  velocity  component  for  approximation  are  a  locally 
specified  primary-flow  direction  and  a  decomposition  of  the  secondary  velocity 
field.  The  local  flow  field  velocity  vector  is  expressed  as  a  primary  flow 
component  in  the  primary-flow  direction  and  a  secondary  flow  component  normal 
to  this  direction.  The  secondary  velocity  component  is  decomposed  into 
components  defined  from  a  scalar-potential  and  a  vector-potential.  The  vector 
potential  defines  the  strearawise  vorticity  and  the  large  secondary  velocity 
component  associated  with  the  streamwise  vorticity.  The  remaining 
scalar-potential  velocity  component  contribution  can  be  assumed  to  be  small  in 
the  transverse  momentum  equations,  but  not  otherwise  neglected.  This  small 
scalar-potential  approximation  along  with  the  viscous  approximation  neglecting 
streamwise  diffusion  is  sufficient  to  establish  a  well-posed  initial  value 
problem,  which  can  be  solved  far  more  economically  than  the  Navier-Stokes 
equations.  No  approximations  are  introduced  for  pressure  in  this  approach. 

Figure  1  illustrates  elements  of  the  velocity  decomposition  and 

approximation.  Figure  la  is  a  schematic  showing  a  local  primary  flow 

direction  and  a  plane  normal  to  the  direction  (local  normal  plane).  The  flow 

field  velocity  vector  is  represented  as  a  primary  flow  (streamwise)  velocity 

component  u  in  the  primary  flow  direction  and  a  secondary  flow  velocity 
-+  P 

vector  V  in  the  local  normal  plane.  The  primary  flow  direction  does  not,  in 
general,  coincide  with  the  marching  coordinate  direction  and  correspondingly 
the  local  normal  plane  does  not  coincide  with  the  transverse  coordinate  plane 
that  is  normal  to  the  marching  coordinate.  The  secondary  flow  velocity 
vector  V  is  decomposed  (Figure  lb)  in  the  local  normal  plane  into  a  component 

-V 

from  a  scalar-potential  (V^ )  and  a  component  from  a  vector  potential 

(V , ).  The  vector  potential  component  V ^ is  associated  with  streamwise 
^  -  >■ 

vorticity  and  Is  expected  to  be  large.  The  scalar  potential  component  is 

the  difference  between  the  total  velocity  component  V  and  and  is  expected 
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to  be  small.  The  small  scalar-potential  approximation  neglects  the  effects  of 
V<j,  in  transverse  momentum  equations.  This  component  is  not  neglected  in  the 
continuity  equation.  A  physical  interpretation  of  the  approximation  for  a 
general  three-dimensional  viscous  flow  field  is  not  easy.  However,  the 
approximate  equations  can  be  reduced  to  two-dimensional  boundary  layer 
and  slender  channel  theory  equations.  In  this  boundary  layer  context,  the 
scalar-potential  component  V ^  corresponds  to  the  small  velocity  component 
normal  to  the  boundary. 

An  important  element  in  the  application  of  the  primary/secondary  flow 
equations  to  compute  three-dimensional  viscous  flow  is  to  identify  a  local 
primary  flow  direction.  Potential  flow  streamlines  conform  to  the  shape  of 
boundaries  and  are  a  suitable  choice  for  the  primary  flow  direction  for  many 
flow  problems.  In  the  case  of  the  submerged  body  flow  field,  the  primary  flow 
direction  was  chosen  to  be  direction  of  potential  flow  streamlines  around  the 
body  at  zero  incidence.  This  potential  flow  was  calculated  using  slender  body 
theory  [5] . 

Validation  of  the  spatial-marching  approach  requires  comparison  of  the 
computation  with  experimental  measurements  and  other  computational  results. 

In  Ref.  fA],  computed  velocities  for  both  laminar  and  turbulent  flow  in  a 
90-degree  square  bend  agree  very  well  with  measurements  of  Taylor,  Whitelaw, 
and  Yianneskis  [6j.  A  comparison  with  a  Navier-Stokes  solution  for  a 
90-degree  channel  bend  gives  excellent  agreement  with  both  pressure  and 
velocity  [A].  In  Ref.  [7],  surface  pressures  computed  for  a  tip  vortex  flow 
field  are  in  good  agreement  with  measurements  of  Gray,  McMahon,  Shenoy  and 
Hammer  [8].  In  the  present  study,  the  flow  past  an  unappended  submarine  hull 
in  drift  has  been  computed  at  a  Reynolds  number  of  16  million  and  incidence  of 
20  degrees. 

The  remainder  of  this  section  contains  a  mathematical  description  of  the 
generalized  primary/secondary  flow  equations  and  the  small  scalar-potential 
approximation.  Details  of  the  derivation  of  these  equations  can  be  found  in 
Ref.  (3).  In  the  summary  that  follows,  i^  denotes  the  unit  vector  basis 
associated  with  the  orthogonal  reference  line  coordinate  system  used  to 
describe  the  equations.  The  vector  ij  is  aligned  with  the  marching 
coordinate  direction  (x  [)  and  vectors  i o  and  i}  Lie  in  the  transverse 
coordinate  plane.  The  vector  basis  e^  and  the  reciprocal  basis  e^  are 
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defined  to  identify  the  local  primary  flow  direction  and  the  local  normal 

plane  for  the  purpose  of  approximating  flow  equations.  The  vector  5^  is 

-2  -3 

aligned  with  the  primary  flow  direction  and  the  vectors  e  and  e  lie  in  the 

“  1c  _ 

local  normal  plane  [Figure  1 . a . ] .  The  vectors  e  and  e,  are  related  and  can 

K 

be  written  in  terms  of  their  components  in  the  coordinate  system  i^  as 

ek*e£  =  ekmem£  =  ^  (2.1) 

The  flow  field  velocity  vector  is  denoted  by  U, 


U  =  u^ij  +  u  2^  2  +  u3*3  =  upe  1  +  ve  2  +  we  3  (2.2) 

The  velocity  component  Up  is  associated  with  the  primary  flow  direction  and 
components  v  and  w  with  the  secondary  flow  field.  The  secondary  flow  velocity 
vector  V  can  be  defined  as 


V  =  ve  2  +  we  3  (2.3) 

The  secondary  flow  velocity  components  are  decomposed  into  scalar  potential 
(v<{,,  w^)  and  vector  potential  components  (v,*,  w  ^) 


V  =  v*  +  V^ 

V  =  ^ 

v  =  v*  a-  v* 

(2.4) 

The  scalar  potential  velocity  components  can  be  defined  in  terms  of  a  scalar 
potential  (  <f>)  and  the  vector  potential  components  in  terms  of  a  vector 
potential  (  <Jj)  . 


vs4>  =  (e^V^ejMfc 


—  ?.  1  - 

Vsx^i!  =  (pe  JV^*e j )i £ 


£  =  2,3  (no  sum) 


(2.5) 


where  P  is  density  and  7g  is  the  surface  operator  defined  as 


7  =  7-7 

s  p 


Vp  =  i^M-9) 


(2.6) 


All  variables  in  the  flow  equations  are  assumed  to  be  nondimensional ,  having 
been  normalized  by  a  suitable  combination  of  reference  quantities.  The 
primary/secondary  flow  equations  are  written  for  the  dependent  variables, 
streaawise  velocity  Up,  scalar  potential  $,  vector  potential  <|>,  vorticity  Q, 
pressure  p,  total  enthalpy  E,  and  density  p. 

Scalar  potential  equation  (continuity): 


7  *(pV  $)  +  7  •  (pe^V*e .  )ii  +  7*(pe^u  )i,  =  0 
s  s  p  j  1  p  k 


(2.7) 


Streaz-ise  momentum  equation: 


pU*7u  -  pU*(U,7)ei  +  ei*Vp  +  ei*Vxy7xU 


(2.8) 


where  Re  is  the  Reynolds  number  defined  from  reference  quantities,  u  the 
normalized  coefficient  of  molecular  viscosity,  and  7  is  an  approximation  to  7 
in  which  second  derivatives  along  the  ij  coordinate  have  been  neglected. 

7ector  potential  equation  (definition  of  vorticity): 


7  x_  (V^i,  )  -  fli,  =0 


(2.9) 
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Vorticlty  transport  equation: 


pup  30 
h  3  xi 


« k  t  —  pur.gp£k___ 

+  7  -p\«  +7  *(e  C,-D  )i  +  (7  -jZ)xlf_  (V-  e,  )i 
s  s  k  l  s  n  3  xj  k  £. 


+  (7  p  V)x [ 7  e£k(V*e,  )]i  +  -i  7  xe  e  -7  xya  i,  =  0 

s  1  s  k  J  £  Re  s  ks  1 


(2.10) 


takes  on 


where  and  D  are  defined  in  (3]  and  V  =  BV^  +  .  The  parameter  8 

values  of  0  or  1,0  =0  identifies  the  small  scalar  potential  approximation. 


Pressure  equation: 


j  k-  -  tk  - 

Vs*(e  ej^VpHjj  +  7  s*  (e  =  0  A  =  2,3  (no  sum)  (2.11) 


is  a  combination  of  terms  from  transverse  momentum  equations  and  is  defined 
in  [3].  The  secondary  flow  velocity  vector,  V,  is  written  as  V  =  BV^  +  in 
F  and  6=0  identifies  the  small  scalar  potential  approximation. 


Energy  equation: 


PU-7E  -  (RePr)"1  7-<7  E  -  (y-1)  M2  Re-1  f -pU-7  xy7  xq 

s  r  1  s  s 


(2.12) 


-  Pr-1  7  kV 
s  s 


qp/2]  =  o 


—  — i  -o  — o 

where  q  =  u  e  +  (6v,  +  v.  )e  +  (Bw,  +  w,  )e  , 

<!>  is  the  dissipation  function,  <  the  coefficient  of  thermal  conductivity,  Pr 
the  Prandtl  number,  and  Mr  the  reference  Mach  number.  Simpler  forms  of  the 
energy  equation  that  can  be  used  at  low  and  moderate  Mach  numbers  are  described 
in  ( 3] . 

State  equat ion  : 


p  =  P  F./y  M^  -  (y  -1  )o  q  2  /2y 


(2.13) 


Equations  (2.7)  to  (2.13),  along  with  the  definition  of  the  velocity 
decomposition  (2.5),  are  the  complete  set  of  primary/secondary  flow  equations. 
The  small  scalar  potential  approximation  (6=0)  makes  these  equations  a 
well-posed  initial  value  problem  along  the  streamwise  coordinate  (x  j). 

3.  SOLUTION  OF  PRIMARY/ SECONDARY  FLOW  EQUATIONS 

The  choice  of  a  solution  algorithm  for  the  primary/secondary  flow 
equations  must  exploit  the  small  scalar  potential  approximation  for  economy. 

The  potential  for  economy  lies  in  the  initial-value  nature  of  the  primary/ 
secondary  flow  equations  and  the  ability  to  solve  such  a  problem  by  a  spatial 
marching  scheme.  Such  a  scheme  will  be  developed  in  this  section  for  the 
primary/secondary  flow  equations. 

An  evolution  equation  can  be  written  in  the  form 

A  —  +  D(  $)  =  0  (3.1) 

3t 

where  <j>  is  a  vector  of  dependent  variables,  t  is  a  time-like  variable  (e.g. 
streamwise  coordinate),  A  is  a  matrix  of  coefficients,  and  D  is  a  multi¬ 
dimensional  differential  operator.  Efficient,  implicit,  noniterative  marching 
schemes  have  been  developed  to  solve  equations  of  the  type  (3.1)  when  the 
matrix  A  is  well-conditioned.  An  example  of  such  a  scheme  is  the  LBI  scheme 
[9]  for  nonlinear  equations  and  is  based  on  the  Douglas-Gunn  ADI  scheme. 
However,  the  matrix  A  for  the  primary/secondary  flow  equations  is  singular  due 
to  the  lack  of  streamwise  derivatives  in  some  of  the  governing  equations.  This 
is  also  true  when  the  same  equations  are  written  in  terms  of  primitive 
variables  [3].  Coupled  implicit  schemes  such  as  the  LBI  scheme  require 
iteration  when  A  is  singular.  This  penalty  in  economy,  due  to  iteration  of  a 
fully  coupled  system  at  each  step,  is  avoided  in  the  present  sequentially 
decoupled  algorithm. 

The  problem  with  the  singularity  of  the  coefficient  matrix  of  the 
streamwise  derivative  is  circumvented  by  decoupling  the  individual  equations  of 
the  system  without  a  streamwise  derivative  and  solving  them  separately  at  each 
streamwise  step  of  a  spatial  marching  scheme  for  the  remaining  equations. 
However,  any  severe  penalty  in  economy  due  to  loss  of  stability  in  the  overall 
scheme  from  the  decoupling  would  be  unacceptable.  The  present  sequentially 
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decoupled  implicit  algorithm  decouples  the  primary/secondary  flow  equations 
into  subsystems  through  sequencing  and  linearization  of  the  subsystems.  Each 
of  the  subsystems  can  then  be  solved  by  economical,  stable,  implicit 
procedures. 

3.1  Sequentially  Decoupled  Implicit  Algorithm 

The  solution  algorithm  for  the  primary/secondary  flow  equations  is  developed  by 
considering  a  more  general  case  of  a  partial  differential  equation  in  the  form 

Mq  =  0  (3.2) 

M  is  a  matrix  operator  (nonlinear)  whose  elements  are  the  coefficients  and 
differential  operators  of  the  individual  equations  of  the  system  and  q  is  a 
vector  of  dependent  variables.  Further,  the  elements  of  q  are  assumed  to  be 
ordered  and  the  sequence  of  equations  in  M  arranged  on  the  basis  of  an 
association  between  an  equation  and  a  dependent  variable.  This  association 
between  individual  equations  of  the  system  and  dependent  variables  will  be 
clarified  with  development  of  the  algorithm.  The  matrix  M  is  decomposed  and 
(3.2)  written  in  the  form 

(L  +  D  +  U)q  =  0  (3.3) 


where  L,  D,  U  are  matrix  operators  containing  the  lower  diagonal,  diagonal  and 
upper  diagonal  elements  of  M,  respectively.  The  matrix  operators  L,  D  and  U 
are  linearized  by  the  following  linearization  scheme 


n,n+l  „  ,  n+1  n. 

-  :  'ij  <’k  ■  \> 


k  =  I ,  .  .  .  (j-l ) 

m=j,  .  .  .N 


(3.  A) 


where  are  the  elements  of  L  and  n  is  the  marching  coordinate  index.  The 
operators  D  and  U  can  be  linearized  in  the  same  manner  as  L.  Equation  (3.3) 
can  be  written  as  a  linear,  decoupled  set  of  equations 


(Dn*n+1  +  Ln’n+1)  qn+1 


+  Un’n+1  q" 


=  0 


(3.3) 
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or 


n,ii+l  n+1  n,n+l  n+1 

D  q  =  -  (L  q 


Un’n+1  q") 


(3.6) 


Equation  (3.6)  is  a  decoupled  system  of  equations  that  can  be  solved  by 
inverting  the  linear  diagonal  operators  in  sequence  using  implicit  schemes. 
Clearly,  the  success  of  the  sequentially  decoupled  algorithm  (3.6)  depends  on 
the  system  of  equations  (3.2)  and  the  choice  of  the  diagonal  matrix  operator, 

D.  The  diagonal  operator,  D,  is  determined  by  the  association  made  between 
individual  equations  of  the  system  and  dependent  variables. 

The  association  of  equations  with  variables  is  straightforward  for  the 
primary/secondary  flow  equation;  ft  with  (2.10),  with  (2.9),  p  with  (2.11), 

Up  with  (2.8),  <J>  with  (2.7),  and  E  with  (2.12).  For  simplicity,  the  implicit 
operator  (Dn,n  *  +  Ln,n  *)  in  (3.5)  for  the  incompressible  form  of  the  primary/ 
secondary  flow  equations  can  be  written  as 


(Dn’n+1  +  Ln’n+1)  qn+1 


-  — 

-  - 

di  i  *  * 

ft 

1  V2 

* 

s 

*31  *32  *  * 

P 

„  „  1  3 

0  0  IT  57 

u 

P 

0  0  0  -1  V2 

3x  s 

* 

_  _ 

_  _ 

1 

where  djj,  d4t4  =  U*V  -  Re-lV  .  The  in  (3.7)  is  used  to  indicate  non-zero 

elements  of  the  explicit  operator,  U  ’  .  The  vector  basis  is  assumed  aligned 

with  coordinate  directions  in  (3.7).  The  vector  basis  introduces  more  complex 
differential  operators  in  individual  equations  but  the  implicit  operator 
retains  essentially  the  same  form  as  (3.7).  Although  no  stability  analysis  of 
(3.6)  has  been  performed,  stability  of  (3.6)  for  the  pr imary/secondary  flow 
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equations  has  been  established  through  extensive  numerical  computation  of 
numerous  test  cases.  Further,  although  the  vorticity  and  vector  potential 
equations  (first  two  equations  in  (3.7))  can  be  uncoupled  in  (3.7),  they  are 
strongly  coupled  at  a  no-slip  boundary  by  boundary  conditions.  Consequently, 
the  vorticity  and  vector  potential  equations  are  solved  as  a  coupled  subsystem 
of  equations  to  avoid  stability  problems.  The  following  sequentially  decoupled 
subsystem  equations  are  obtained  for  the  compressible  primary/secondary  flow 
equations 

1)  vorticity  (2.10),  vector  potential  (2.9)  for  ft,  \jj 

2)  pressure  equation  (2.11)  for  p 

3)  energy  equation  (2.12)  for  E 

4)  streamwise  momentum  (2.8)  for  Up 

3)  equation  of  state  (2.13)  for  p  (algebraic) 

6)  scalar  potential  equation  (2.7)  for  <j> 

The  sequence  completes  the  solution  of  the  primary/secondary  flow  equations  for 
one  streamwise  step  (n+1)  from  the  solution  at  the  previously  computed  step  (n) 
(initial  conditions  for  n=l).  Each  of  the  individual  subsystem  equations  can 
be  solved  by  a  convenient,  efficient  procedure.  Numerical  procedures  used  in 
the  present  work  for  solving  decoupled  subsystem  equations  are  described  in  a 
later  section. 

3.2  Differencing  Procedures  and  Boundary  Conditions 

The  governing  equations  (2.7) — ( 2.12)  are  transformed  to  a  body-fitted 
coordinate  system  in  which  the  equations  are  discretized.  Body-fitted 
coordinate  transformations  suitable  for  the  primary/secondary  flow  equations 
are  discussed  in  (3).  The  discretization  employs  two-point  backward 
differences  for  streamwise  derivatives  and  three-point  central  differences  for 
trans\erse  derivatives.  The  body-fitted  coordinate  transformation  allows  a 
variable  streamwise  step  size  and  concentration  of  grid  points  in  regions  of 
expected  large  gradients.  A  centered  Crank-Nico Ison  formulation  could  be  used 
for  streamwise  derivatives  but  the  formal  accuracy  of  the  scheme  would  not  be 
improved  because  of  the  first-order  accurate  linearization  scheme  (3.4)  used  to 
decouple  subsystem  equations. 
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Implicit  boundary  conditions  are  used  in  all  subsystem  equations. 

Boundary  conditions  that  are  not  straightforward  are  no-slip  conditions  on  a 

solid  boundary  and  boundary  conditions  on  pressure.  These  conditions  are 

described  here.  The  streamwise  velocity  Up  is  specified  as  zero  in  the 

streamwise  momentum  equation  (2.8)  on  a  no-slip  boundary.  The  condition 

v  =  w  =  0  must  be  expressed  in  terms  of  $ ,  t|»  and  Q.  The  normal  component  of 

the  scalar  potential  velocity  is  specified  as  zero  by  i  *V  =  0,  where  i  is 

n  $  n 

the  unit  normal  vector  to  the  boundary.  V  is  specified  in  terms  of  <$>  to 

9 

obtain  a  Neumann  condition  for  the  scalar  potential  equation  (2.7).  The  normal 
component  of  the  vector  potential  velocity  is  specified  zero  by  prescribing  ^  = 
0  on  the  boundary.  The  tangential  component  of  the  transverse  velocity  is 
specified  at  the  no-slip  boundary  by 

v[(v* +  V*2  +  (w* +  =  0  (3-8) 

where  it  is  the  unit  tangent  vector  to  the  boundary.  Equation  (3.8)  combined 
with  equation  (2.9)  is  used  to  obtain  ft  on  the  boundary  in  terms  of  4>  and  <J>. 
This  relation  for  ft,  along  with  prescribing  ip  =  0,  provides  coupled  boundary 
conditions  for  the  vector  potential  (2.9)  and  vorticity  equations  (2.10).  The 
coupled  boundary  conditions  require  the  equations  to  be  solved  as  a  coupled 
subsystem  from  stability  considerations. 

Boundary  conditions  for  pressure  are  obtained  by  computing  the  normal 
gradient  of  pressure  on  the  boundary  from  transverse  momentum  equations.  The 
prescription  of  the  normal  gradient  of  pressure  in  the  pressure  equation  (2.11) 
results  in  a  Neumann  problem.  Solution  of  the  Neumann  problem  leaves  an 
undetermined  constant  in  the  transverse  plane  that  could,  however,  be  a 
function  of  the  streamwise  coordinate.  Pressure,  p,  is  written  as 

p  =  Pm(xi)  +  Py  (*i i  x2 ,  x3)  (3.9) 

where  pm  is  the  undetermined  constant  of  the  Neumann  problem  and  pv  is  the 
solution  of  the  pressure  equation.  The  value  of  pm  is  required  in  the 
solution  of  the  streamwise  momentum  equation  (2.8)  in  the  form  of  the 
streamwise  gradient,  .  In  external  flows,  pro  can  be  determined  by  equating 
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the  right-hand  side  of  (3.9)  to  prescribed  pressure  at  a  point  in  freestream  in 
each  transverse  plane.  Such  a  condition  is  not  available  in  internal  flows  and 
pm  has  to  be  determined  to  ensure  that  the  integral  mass  flux  relation 


{  pU-ijdA  =  Q  =  constant  (3.10) 

A 


is  satisfied.  For  internal  flows,  pm  is  adjusted  by  a  secant  iteration 
during  solution  of  the  streamwise  momentum  equation  (2.8)  so  that  integral  mass 
flux  relation  (3.10)  is  satisfied. 

Far  field  boundary  conditions  for  velocity  for  the  submerged  body  flow 
field  are  specified  in  terms  of  freestream  conditions  and  the  angle  of  attack. 
The  streamwise  velocity  and  vorticity  are  extrapolated  at  the  outer  boundary. 
The  vector  potential  component  of  the  transverse  velocity  field  is  specified  in 
terms  of  the  angle  of  attack  (o )  as 

i2*Vg3fl|»ii  =  U^sino  (3.11) 

Equation  (3.11)  is  integrated  along  the  outer  boundary  for  a  specification  of 
i(>.  The  scalar  potential  component  of  the  transverse  velocity  is  specified  in 
terras  of  the  normal  gradient  of  the  scalar  potential 

Q.  -  r.  - 

i  «V  <f>  =  ~l  _J-i«{£2-  +  _i:  +  -5°)  <5x26x3  (3.12) 

ns  .  2ir  n  r  r .  r 

3  1  0 

where  r  is  the  vector  from  the  grid  point  j  to  the  boundary  point  under 
consideration,  is  the  vector  from  the  image  of  the  grid  point  j  to  the 
boundary  point,  and  tq  is  the  vector  from  the  origin  to  the  boundary  point. 

Qj  is  the  source  strength  obtained  from  (2.7)  as 

Q  =7  -(oe'Ve  )ii  +V*(pe'C1u  )i,  (3.13) 

J  P  J  P  k 

The  summation  in  (3.12)  is  carried  out  over  all  grid  points  in  the  transverse 
plane,  extending  2?  around  the  body. 
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3.3  Turbulence  Model 


The  mixing  length  turbulence  model  developed  by  Baldwin  and  Lomax  [10]  was 
utilized  to  compute  the  eddy  viscosity.  In  utilizing  the  model  for  the 
submerged  body  flow  field,  problems  were  encountered  in  the  computation  of  the 
vorticity  function  (FMAX).  In  the  region  of  the  lee-side  vortex,  a  second  peak 
in  the  vorticity  function  was  found  near  the  center  of  the  vortex.  This  peak 
was  of  greater  value  than  the  first  peak  found  in  the  boundary  layer  region  of 
the  flow  field.  If  the  second  peak  was  utilized  in  the  model,  large  eddy 
viscosities  were  computed  in  the  region  of  the  lee-side  vortex  which  is 
physically  a  region  of  rotational  inviscid  flow.  Large  eddy  viscosities  in  the 
region  of  the  vortex  dissipated  the  vortex  in  the  computations.  To  overcome 
the  problem,  the  turbulence  model  was  modified  to  restrict  the  search  for  the 
peak  in  the  vorticity  function  to  y+  =  2000.  Similar  modifications  to  the 
mixing  length  model  have  been  used  in  the  computations  of  Degani  and 
Schiff  [11]  and  Panaras  and  Steger  [12]. 

3.4  Summary 

A  summary  of  the  procedure  used  to  advance  the  solution  a  single 
streanwise  step  to  the  (n+1)  level  from  previously  computed  quantities  at  the 
n  level  follows.  Coordinate  and  vector  basis  information  required  to  solve  the 
primary/secondary  equations  are  prescribed  input.  Governing  equations  are 
assumed  to  be  linearized  and  discretized  as  previously  described. 

(1)  Equations  (2.9)  and  (2.10)  form  a  linear  coupled  system  for  i|0+l 
and  fin+1 ,  which  is  solved  as  a  (2*2)  coupled  system.  For  this 
purpose,  artificial  time  derivatives  are  added  to  each  equation  and 
an  iterative  block-implicit  scheme  [9]  is  used.  In  prescribing 
no-slip  boundary  conditions,  the  tangential  component  (3.8)  contains 
a  contribution  from  4>;  this  contribution  is  evaluated  using  $8. 
Linearization  (3.4)  provides  that  terms  involving  other  unknown 
dependent  variables  are  evaluated  at  the  n-level. 
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(2)  The  pressure  equation  (2.11)  is  solved  for  p|J+^  using  an  iterative 
scalar  ADI  scheme.  Updated  values  of  vI)+*  and  w7+^  are  available 

<P  <l> 

while  other  dependent  variables  are  evaluated  at  the  n-level. 

(3)  The  energy  equation  (2.12)  is  solved  for  En_t  ^  using  a  scalar  ADI 
scheme. 

(4a)  The  streamwise  momentum  equation  (2.8)  is  solved  for  ujj+^  using  a 

scalar  ADI  scheme.  The  streamwise  pressure  gradient  is  evaluated  at 
the  (n+1)  level. 

(4b)  The  density  pn+1  is  evaluated  algebraically  from  the  equation  of 
state  (2.13) 

(5)  Finally,  the  scalar  potential  equation  (2.7)  is  solved  for  <fin+l 
using  an  iterative  scalar  ADI  scheme. 

4.  COMPUTED  RESULTS  FOR  FLOW  PAST  A  SUBMERGED  BODY  IN  DRIFT 

The  primary  focus  of  the  present  effort  has  been  to  demonstrate  the 
capability  to  compute  the  three-dimensional  viscous  flow  around  an  unappended 
submarine  hull  in  drift,  utilizing  the  generalized  primary/secondary  flow 
equations.  Results  of  the  computed  flow  field  around  an  unappended  submarine 
hull  in  drift  (20  degrees)  are  presented  here.  These  results  are  compared  with 
the  limited  experimental  data  available  for  this  flow  case.  The  Reynolds 
number  of  the  flow  based  on  the  length  of  the  hull  was  1.68 *10  7.  Computations 
were  carried  out  for  one-half  of  the  symmetric  flow  field  on  a  grid  of  99 
circumferential  points,  70  radial  points  and  191  streamwise  points  (1.32 *10& 
total  grid  points). 

Figure  2  provides  a  visualization  of  the  roll-up  of  the  lee-side  vortex 
through  particle  paths  traced  in  the  computed  flow  field.  Particle  traces  are 
started  in  the  upstream  boundary  layer  at  four  streamwise  stations.  Figure  3 
shows  a  contour  plot  of  the  streamwise  velocity  at  several  transverse  planes. 
Figure  4  shows  a  similar  contour  plot  of  the  computed  static  pressure 
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at  several  transverse  planes.  Figures  2-4  provide  a  basis  for  the  description 
of  the  physical  processes  in  the  flow  around  the  hull  in  drift.  The  cross  flow 
component  of  the  velocity  field,  due  to  drift  of  the  hull,  convects  low 
momentum  fluid  in  the  boundary  layer  on  the  hull  towards  the  lee  side.  This  is 
clearly  visible  in  the  trajectories  of  the  particle  paths  in  the  upstream  flow 

(Figure  2).  The  low  momentum  fluid  on  the  lee  side  rolls  up  into  a  vortex  and 

is  convected  downstream  by  the  streamwise  velocity.  The  upstream  static 
pressure  distribution  on  the  hull  (Figure  4)  shows  a  low  pressure  region  (blue) 
on  top  of  the  hull  caused  by  the  acceleration  of  the  cross  flow  velocity 
component.  The  pressure  recovers  on  the  lee  side  of  the  body.  However,  as  the 
lee  side  vortex  is  generated,  a  low  pressure  region  forms  at  the  center  of  the 

vortex.  The  low  pressure  region  in  the  vortex  core  significantly  distorts  the 

pressure  field  at  the  stern.  The  high  pressure  (red)  region  at  the  nose  of  the 
body  is  due  to  the  stagnation  point  at  the  nose.  The  stagnation  point  is  also 
the  source  of  the  streamwise  velocity  gradients  at  the  nose  (Figure  3). 

Figures  2  and  3  also  show  boundary  layer  fluid  being  convected  away  from  the 
nose,  on  the  lee  side,  along  the  symmetry  plane.  It  is  not  yet  clear  whether 
this  aspect  of  the  flow  field  has  been  observed  in  experiments  or  is  related  to 
the  initial  velocity  profiles  assumed  in  the  computations. 

Figures  5  and  6  show  vector  plots  of  the  transverse  velocity  field  at  two 
streamwise  stations  (x  =  0.8  and  0.95).  These  plots  provide  a  visualization  of 
the  structure  of  the  lee  side  vortex.  Transverse  velocity  vectors  have  been 
colored  by  the  static  pressure  field,  providing  a  relation  between  the 
transverse  velocity  and  pressure  fields.  The  plots  clearly  show  the  low 
pressure  regions  at  the  vortex  center  and  on  the  hull.  The  pressure  field  also 
shows  the  interaction  between  the  pressure  field  due  to  the  vortex  and  the  body 
pressure  field.  Figure  7  compares  the  computed  transverse  velocity  field  at 
x  =  0.95  with  measurements  at  x  =  1.0.  These  measurements  were  provided  by 
DTRC.  The  strength  of  the  computed  and  measured  vortices  compare  well  with 
each  other  but  the  computed  vortex  is  located  closer  to  the  body  than  in  the 
measurements.  The  source  of  discrepancy  in  the  location  of  the  vortex  could  be 
due  to  modifications  adopted  in  the  turbulence  model  to  compute  the  boundary 
layer  length  scale  or  due  to  inadequate  resolution  of  the  computed  flow  field 
in  the  stern  region.  Some  of  the  discrepancy  could  also  be  attributed  due  to 
the  different  streamwise  locations  of  the  measurements  and  computations.  The 
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computations  were  not  carried  out  to  x=1.0  due  to  the  computational  grid 


resolution  being  inadequate  in  this  region,  in  the  present  case.  Further,  the 
measured  lee  side  vortices  are  asymmetric.  Means  to  include  this  effect  in  the 
computations  were  beyond  the  scope  of  the  present  effort. 

Figure  8  shows  the  distribution  of  the  computed  axial  and  side  forces  on 
the  hull.  The  figure  shows  the  distribution  of  forces  due  to  the  pressure 
field  alone  and  the  distribution  due  to  the  total  force  (pressure  and  shear). 
The  viscous  shear  stress  does  not  alter  the  distribution  of  the  force  but  does 
change  the  magnitude  of  the  force  slightly.  As  would  be  expected,  the  viscous 
shear  stress  affects  the  axial  force  to  a  greater  magnitude  than  the  side 
forces.  Although  the  viscous  shear  stress  does  not  significantly  affect  the 
side  forces,  the  pressure  field  in  itself  is  significantly  altered  by  the 
viscous  flow,  as  compared  to  potential  flow,  producing  a  net  side  force  and 
axial  force  on  the  hull.  The  side  force  distribution  on  the  hull  (Figure  8) 
shows  a  positive  force  distribution  at  the  nose  which  decreases  beyond  the 
nose.  The  side  force  increases  again  as  the  lee  side  vortex  is  generated.  The 
side  force  is  negative  at  the  stern.  The  side  force  distribution  produces  a 
net  positive  force  on  the  hull.  The  axial  force  distribution  shows  small 
regions  of  negative  drag  induced  by  the  pressure  field  at  the  nose  and  stern  of 
the  hull.  The  axial  force  distribution  produces  a  net  positive  drag  on  the 
hull,  as  would  be  expected. 


5.  CONCLUSIONS 

The  feasibility  of  computing  the  three-dimensional  high  Reynolds  number 
flow  around  unappended  hulls  in  drift  using  the  generalized  primary/secondary 
flow  equations  and  the  spatial  marching  computation  procedure  has  been 
demonstrated.  Important  features  of  the  flow  field  have  been  captured  in  the 
computations.  The  computations  compare  well,  in  terms  of  the  vortex  strength, 
with  the  limited  experimental  data  available  for  the  flow  case.  Computed 
forces  on  the  hull  show  expected  features.  The  computation  procedure  is  very 
economical  considering  details  of  the  flow  field  provided.  The  procedure 
typically  requires  20  minutes  of  CRAY  X-MP  CPU  time  for  10^  grid  points. 
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6.  FUTURE  WORK 


Computed  flow  field  results  need  more  detailed  validation  by  comparison 
with  experimental  data  to  identify  any  sources  of  discrepancy  in  the 
computations.  Determination  of  the  length  scale  for  use  in  the  mixing  length 
turbulence  model  for  the  three-dimensional  viscous  flow  around  the  hull  needs 
to  be  carefully  evaluated.  A  two-equation  turbulence  model  may  be  required  in 
the  computations.  The  computation  procedure  needs  to  be  extended  to  include 
the  effects  of  appendages  on  the  hull.  The  spatial  marching  procedure  could 
also  be  extended  to  include  unsteady  maneuvers.  Such  a  procedure  would  still 
provide  significant  computational  advantages  over  solution  of  the  unsteady 
Navier-Stokes  equations,  especially  for  the  long  time  scales  involved  in  a 
maneuver  and  for  low  mach  number  flows.  The  present  procedure  could  also  be 
used  as  a  tool  to  provide  insight  into  the  basic  physical  processes  of 
asymmetric  flows  around  the  hull. 
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Figure  lb.  Vector  Decomposition  of  Transverse  Velocity. 
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Figure  Captions  for  Color  Plates 


Figure  2. 
Figure  3. 

Figure  A. 
Figure  5. 

Figure  6. 


Particle  paths  showing  lee-side  vortex  roll-up. 

Contour  plot  of  streamwise  velocity  at  several  transverse 
planes . 

Contour  plot  of  static  pressure  at  several  transverse  planes. 

Vector  plot  of  the  transverse  velocity  field  colored  by  static 
pressure  (x  =  0.8). 

Vector  plot  of  the  transverse  velocity  field  colored  by  static 
pressure  (x  =  0.95)  . 
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Figure  7.  Comparison  of  Computed  and  Measured  Transverse  Velocity  Field. 


